package test.demo;

import com.alibaba.fastjson.JSON;
import tech.waterism.modelbase.PredictBase;
import tech.waterism.modelbase.PredictResult;

import java.util.*;

/**
 * 新安江模型（新安江产汇流）
 *
 * @author 初京刚
 * @version 1.0
 * @date 2022-1-7
 */
public class XAJ extends PredictBase {

    //region Initial conditions, Parameters, Inputs, Statuses and Outputs

    //region Initial conditions for dunneMechanism

    /**
     * 流域上层、下层、深层张力水蓄量初值、初始自由水蓄量、初始产流面积
     */
    private double[] WInitialDay;

    /**
     * 流域上层张力水含量
     */
    private double wup;

    /**
     * 下层张力水含量
     */
    private double wlp;

    /**
     * 深层张力水含量
     */
    private double wdp;

    //endregion

    //region Initial conditions for freeTank

    /**
     * 自由水水深
     */
    private double sp;

    /**
     * 产流面积系数
     */
    private double frp;

    //endregion

    //region Initial conditions for LinerReservoir

    /**
     * 地表径流流量、壤中流流量、地下径流流量
     */
    private double[] QInitial;

    /**
     * 地表径流流量
     */
    private double qsp;

    /**
     * 壤中流流量
     */
    private double qip;

    /**
     * 地下径流流量
     */
    private double qgp;

    //endregion

    //region Initial conditions for TimeLag

    /**
     * 初始总基流
     */
    private double[] QSInitial;

    //endregion

    //region Initial conditions for muskingum

    /**
     * 马法初始流量，河段下断面初始流量
     */
    private double[] QXInitial;

    //endregion

    //region Parameters for dunneMechanism

    /**
     * 流域平均张力水容量（mm）、流域蓄水容量
     */
    private double Wm;

    /**
     * 上层张力水容量（mm）
     */
    private double Um;

    /**
     * 下层张力水容量（mm）
     */
    private double Lm;

    /**
     * 深层张力水容量
     */
    private double Dm;

    /**
     * 流域蒸散发折算系数
     */
    private double Kc;

    /**
     * 张力水蓄水容量曲线的方次、蓄水容量曲线指数
     */
    private double B;

    /**
     * 深层蒸散发折算系数
     */
    private double C;

    /**
     * 不透水的面积占全流域面积的比
     */
    private double Im;

    //endregion

    //region Parameters for freeTank

    /**
     * 流域自由水蓄水容量
     */
    private double SM;

    /**
     * 自由水蓄水容量曲线指数
     */
    private double EX;

    /**
     * 地下水径流出流系数（cjg）表层自由水蓄水库对地下水的日出流系数--有Kg+Ki=0.7
     */
    private double KG;

    /**
     * 壤中水径流出流系数
     */
    private double KI;

    //endregion

    //region Parameters for LinerReservoir

    /**
     * 壤中流消退系数
     */
    private double CI;

    /**
     * 地下径流消退系数
     */
    private double CG;

    //endregion

    //region Parameters for TimeLag

    /**
     * 地面径流消退系数
     */
    private double CS;

    /**
     * 滞时（h）
     */
    private double L;

    //endregion

    //region Parameters for muskingum

    /**
     * 子河段洪水波传播时间
     */
    private double KE;

    /**
     * 河段数
     */
    private int N;

    /**
     * 子河段流量比重因子
     */
    private double XE;

    //endregion

    //region Inputs for dunneMechanism

    /**
     * 流域 月蒸发量（mm）
     */
    private double[] ES;

    /**
     * 流域 时段蒸发量（mm）
     */
    private double[] floodE;

    /**
     * 降水时间
     */
    private Date[] floodTm;

    /**
     * 降水量
     */
    private double[] floodDrp;

    /**
     * 降水序列长度
     */
    private Integer floodRainRange;

    /**
     * 降水时段长
     */
    private double timeInterval;

    //endregion

    //region Inputs for LinerReservoir

    /**
     * 流域面积
     */
    private double Area;

    //endregion

    //region Statuses for dunneMechanism

    /**
     * 上层张力水含量
     */
    private double[] Wu;

    /**
     * 下层张力水含量
     */
    private double[] Wl;

    /**
     * 深层张力水含量
     */
    private double[] Wd;

    //endregion

    //region Statuses for freeTank

    /**
     * 产流面积自由水水深
     */
    private double[] S;

    /**
     * 产流面积系数
     */
    private double[] FR;

    //endregion

    //region Statuses for LinerReservoir

    /**
     * 坡面汇流-地表径流状态
     */
    private double[] QSS0;

    /**
     * 坡面汇流-壤中流状态
     */
    private double[] QIS0;

    /**
     * 坡面汇流-地下径流状态
     */
    private double[] QGS0;

    //endregion

    //region Statuses for TimeLag

    /**
     * 初始总基流
     */
    private double[][] QS0;

    //endregion

    //region Statuses for muskingum

    /**
     * 马法初始流量，河段下断面各时段初始流量
     */
    private double[][] QX0;

    //endregion

    //region Outputs for dunneMechanism

    /**
     * 总产流
     */
    private double[] R;

    /**
     * 净雨
     */
    private double[] PE;

    /**
     * 不透水面积产流（测试）
     */
    private double[] RIm;

    //endregion

    //region Outputs for freeTank

    /**
     * 分水源-地表水
     */
    private double[] RS;

    /**
     * 分水源-壤中水
     */
    private double[] RI;

    /**
     * 分水源-地下水
     */
    private double[] RG;

    /**
     * 模拟分水源总产流
     */
    private double[] RunoffSim;

    //endregion

    //region Outputs for LinerReservoir

    /**
     * 坡面汇流-地表径流
     */
    private double[] QS;

    /**
     * 坡面汇流-壤中流
     */
    private double[] QI;

    /**
     * 坡面汇流-地下径流
     */
    private double[] QG;

    //endregion

    //endregion

    public XAJ(Map<String, String> data) {
        super();
        wup = Double.parseDouble(data.get("WUP"));
        wlp = Double.parseDouble(data.get("WLP"));
        wdp = Double.parseDouble(data.get("WDP"));
        sp = Double.parseDouble(data.get("SP"));
        frp = Double.parseDouble(data.get("FRP"));
        WInitialDay = new double[]{wup, wlp, wdp, sp, frp};
        QSInitial = JSON.parseObject(data.get("QSIG"), double[].class);
        QXInitial = JSON.parseObject(data.get("QXSIG"), double[].class);
        qsp = Double.parseDouble(data.get("QSP"));
        qip = Double.parseDouble(data.get("QIP"));
        qgp = Double.parseDouble(data.get("QGP"));
        QInitial = new double[]{qsp, qip, qgp};
        Wm = Double.parseDouble(data.get("WM"));
        Um = Double.parseDouble(data.get("WUMx")) * Wm;
        Lm = Double.parseDouble(data.get("WLMx")) * (Wm - Um);
        Dm = Wm - Um - Lm;
        Kc = Double.parseDouble(data.get("K"));
        B = Double.parseDouble(data.get("B"));
        C = Double.parseDouble(data.get("C"));
        Im = Double.parseDouble(data.get("IM"));
        SM = Double.parseDouble(data.get("SM"));
        EX = Double.parseDouble(data.get("EX"));
        KG = Double.parseDouble(data.get("KG"));
        KI = Double.parseDouble(data.get("KI"));
        ES = JSON.parseObject(data.get("ES"), double[].class);
        floodTm = JSON.parseObject(data.get("dt"), Date[].class);// input[i].data.dt
        floodDrp = JSON.parseObject(data.get("rain"), double[].class);// input[i].data.rain
        floodE = JSON.parseObject(data.get("evaporation"), double[].class);  // input[i].data.evaporation
        floodRainRange = floodDrp.length;
        timeInterval = Double.parseDouble(data.get("clen"));                  // index.clen
        RunoffSim = new double[floodRainRange];
        R = new double[floodRainRange];
        PE = new double[floodRainRange];
        RIm = new double[floodRainRange];
        RS = new double[floodRainRange];
        RI = new double[floodRainRange];
        RG = new double[floodRainRange];
        QSS0 = new double[floodRainRange + 1];
        QIS0 = new double[floodRainRange + 1];
        QGS0 = new double[floodRainRange + 1];
        Area = Double.parseDouble(data.get("F"));
        CS = Double.parseDouble(data.get("CS"));
        CI = Double.parseDouble(data.get("CI"));
        CG = Double.parseDouble(data.get("CG"));
        L = Double.parseDouble(data.get("LAG"));
        KE = Double.parseDouble(data.get("KK"));
        XE = Double.parseDouble(data.get("X"));
        N = Integer.parseInt(data.get("MP"));

        QS = new double[floodRainRange];
        QI = new double[floodRainRange];
        QG = new double[floodRainRange];
    }

    /**
     * 模型计算
     */
    public Map<String, double[]> XAJmodel() {

        //蓄满产流
        dunneMechanism();

        // 分水源
        freeTank();

        for (Integer j = 0; j < floodRainRange; j++) {
            RunoffSim[j] = RS[j] + RI[j] + RG[j];
        }

        // 坡面汇流-线性水库
        LinerReservoir();

        double[] Qtmp = new double[floodRainRange];//单元面积河网总入流（m3/s）

        for (Integer j = 0; j < floodRainRange; j++) {
            Qtmp[j] = QS[j] + QI[j] + QG[j];
            if (Qtmp[j] < 0.0) {
                Qtmp[j] = 0.0;
            }
        }

        //河网汇流-滞后演算法,单元面积出口流量（m3/s）
        double[] Qf = TimeLag(Qtmp);

        // 河道演进
        double[] QSim = muskingum(Qf);

        Map<String, double[]> simulationResult = new HashMap<>();
        simulationResult.put("rain", floodDrp);
        simulationResult.put("PE", PE);
        simulationResult.put("runoffSim", RunoffSim);
        simulationResult.put("RS", RS);
        simulationResult.put("RI", RI);
        simulationResult.put("RG", RG);
        simulationResult.put("QS", QS);
        simulationResult.put("QI", QI);
        simulationResult.put("QG", QG);
        simulationResult.put("QSim", QSim);
        return simulationResult;
    }

    //    /**
    //     * 对降雨量的隶属度计算
    //     *
    //     * @param rain 降水量
    //     * @return 降雨对小雨、中雨、大雨的隶属度
    //     */
    //    public double[] calcRainMembership(double rain) {
    //        double[] Pms = {5.0, 17.5, 30.0};
    //
    //        double[] Ms = new double[Pms.length]; // 该时段降雨对小、中、大雨的隶属度
    //
    //        if (rain <= Pms[0]) {
    //            Ms[0] = 1.0; // 对小雨的隶属度
    //            Ms[1] = 0.0; // 对中雨的隶属度
    //            Ms[2] = 0.0; // 对大雨的隶属度
    //        } else if (rain > Pms[0] && rain <= Pms[1]) {
    //            Ms[0] = (Pms[1] - rain) / (Pms[1] - Pms[0]);
    //            Ms[1] = (rain - Pms[0]) / (Pms[1] - Pms[0]);
    //            Ms[2] = 0.0;
    //        } else if (rain > Pms[1] && rain <= Pms[2]) {
    //            Ms[0] = 0.0;
    //            Ms[1] = (Pms[2] - rain) / (Pms[2] - Pms[1]);
    //            Ms[2] = (rain - Pms[1]) / (Pms[2] - Pms[1]);
    //        } else {
    //            Ms[0] = 0.0;
    //            Ms[1] = 0.0;
    //            Ms[2] = 1.0;
    //        }
    //
    //        return Ms;
    //    }

    /**
     * 蒸发计算
     *
     * @param date   降水时间
     * @param timeIn 时间间隔 (h)
     * @return 蒸发
     */
    private double calcEvaporationPotential(Date date, double timeIn) {
        Calendar time = Calendar.getInstance();
        time.setTime(date);
        Integer Month = time.get(Calendar.MONTH) + 1; // 获得降雨开始月份

        /*
        double[] Ev = new double[3]; // 根据月份确定该时刻小雨、中雨、大雨的蒸发能力;

        if (Month == 4)
            Ev = Eall[0];
        else if (Month == 5)
            Ev = Eall[1];
        else if (Month == 6)
            Ev = Eall[2];
        else if (Month == 7)
            Ev = Eall[3];
        else if (Month == 8)
            Ev = Eall[4];
        else if (Month == 9)
            Ev = Eall[5];
        else if (Month == 10)
            Ev = Eall[6];
        else
            System.out.println("The Evaporation of some MONTHs is not found.");

        rain = rain * 24 / timeIn;  //时段降雨量转换成日降雨量

        double[] Ms = calcRainMembership(rain); // 计算该场降雨对小雨、中雨、大雨的隶属度

        */

        double Ep = ES[Month - 1] / 30; // 月蒸发量转换成日蒸发量

        /*
        for (int i = 0; i < Ms.length; i++) { // 隶属度和蒸发数组加权求和
            Ep += Ev[i] * Ms[i];
        }
        */

        Ep = Ep * timeIn / 24; // 日蒸发量转换成时段蒸发量
        return Ep;
    }

    /**
     * 蓄满产流模型
     *
     * @return 总产流、净雨
     */
    private void dunneMechanism() {
        double Wmm = Wm * (1.0 + B) / (1.0 - Im); // 流域单点最大蓄水容量 <Computer Models of Watershed Hydrology> 219页公式7.4、7.5
        Wu = new double[floodRainRange + 1]; // 上层张力水蓄量，随时段变化
        Wl = new double[floodRainRange + 1]; // 下层张力水蓄量，随时段变化
        Wd = new double[floodRainRange + 1]; // 深层张力水蓄量，随时段变化

        Wu[0] = WInitialDay[0]; // 上层张力水蓄量初值
        Wl[0] = WInitialDay[1]; // 下层张力水蓄量初值
        Wd[0] = WInitialDay[2]; // 深层张力水蓄量初值

        double[] Eu = new double[floodRainRange]; // 上层蒸散发量，随时段变化
        double[] El = new double[floodRainRange]; // 下层蒸散发量，随时段变化
        double[] Ed = new double[floodRainRange]; // 深层蒸散发量，随时段变化

        double EPt; // EPt为蒸散发，等于Kc * E（实测蒸发，即ES[]或时段）
        double Wmt; // 时段初始的流域平均张力水蓄量
        double At; // 与Wmt值相对应的纵坐标值（流域单点最大蓄水量）
        double WPt; //时段末上层张力水蓄量
        double WWPt;//时段末上层和下层张力水蓄量;

        if (Dm < 0) {
            System.out.println("Dm:" + Dm);
            Dm = 0.0;
        }
        for (Integer i = 0; i < floodRainRange; i++) {
            Wmt = Wu[i] + Wl[i] + Wd[i];
            if (Wmt < 0.0) {
                System.out.println("W:" + Wmt);
            }
            if (Wmt > Wm) {
                Wmt = Wm;
            }
            //蒸发计算
            if (floodE == null || floodE.length != floodRainRange) {
                EPt = Kc * calcEvaporationPotential(floodTm[i], timeInterval);
            } else {
                EPt = Kc * floodE[i];
            }
            //EPt = Kc * calcEvaporationPotential(floodTm[i], floodDrp[i], timeInterval);

            // PE[i]为降雨减去蒸发，可正可负
            PE[i] = floodDrp[i] - EPt;

            //如果有净雨产生，上、下、深层张力水蓄量主要受到降雨影响，是加法 *********/
            if (PE[i] >= 0.0) {
                //计算总产流
                At = Wmm * (1.0 - Math.pow((1.0 - Wmt / Wm), 1.0 / (1.0 + B)));
                if (PE[i] + At <= Wmm) {
                    R[i] = PE[i] + Wmt - Wm + Wm * Math.pow((1 - (PE[i] + At) / Wmm), 1 + B);
                } else {
                    R[i] = PE[i] + Wmt - Wm;
                }

                // 不透水面积上直接产流
                RIm[i] = PE[i] * Im;

                //各层蒸散发初值
                Eu[i] = EPt;
                El[i] = 0.0;
                Ed[i] = 0.0;

                /********* 分情况讨论上、下、深层张力水蓄量是否饱和，计算时段末Wu、Wl、Wd *********/
                // 如果上层未饱和
                // PE[i]里面已经考虑蒸发了
                WPt = Wu[i] + PE[i] - R[i];
                if (WPt < Um) {

                    /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                    //                        if (i + 1 < floodRainRange) {
                    Wu[i + 1] = WPt;
                    Wl[i + 1] = Wl[i];
                    Wd[i + 1] = Wd[i];
                    //                        }
                }
                // 如果上层饱和，分情况讨论下层是否饱和
                else {
                    // 在上层饱和前提下，如果下层未饱和
                    WWPt = Wu[i] + Wl[i] + PE[i] - R[i];

                    /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                    //                        if (i + 1 < floodRainRange) {
                    if (WWPt < Um + Lm) {
                        Wu[i + 1] = Um;
                        Wl[i + 1] = WWPt - Um;
                        Wd[i + 1] = Wd[i];
                    }
                    // 在上层饱和前提下，如果下层饱和
                    else {
                        Wu[i + 1] = Um;
                        Wl[i + 1] = Lm;
                        Wd[i + 1] = Wmt + PE[i] - R[i] - Wu[i + 1] - Wl[i + 1];
                    }
                    //                        }
                }
            }
            //如果没净雨产生，上、下、深层张力水蓄量主要受到蒸发影响，是减法 *********/
            else if (PE[i] < 0.0) {

                //总产流为0
                R[i] = 0.0;

                // 不透水面积上直接产流为0
                RIm[i] = 0.0;

                /********* 分情况讨论上、下、深层蒸散发是否达到能力，计算时段末Wu、Wl、Wd *********/
                // 如果超过蒸散发能力，即Wu[i]+floodDrp[i]-EPt>=0，仅上层有蒸发
                if (Wu[i] + PE[i] >= 0.0) {
                    Eu[i] = EPt;
                    El[i] = 0.0;
                    Ed[i] = 0.0;

                    /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                    //                        if (i + 1 < floodRainRange) {
                    Wu[i + 1] = Wu[i] + PE[i]; // PE[i]=floodDrp[i]-EPt，这里PE[i]是负值
                    Wl[i + 1] = Wl[i];
                    Wd[i + 1] = Wd[i];
                    //                        }
                }
                // 如果没超过蒸散发能力，即Wu[i]+floodDrp[i]-EPt<0，分情况讨论是否超过下层蒸散发能力
                else {
                    Eu[i] = Wu[i] + floodDrp[i];

                    /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                    //                        if (i + 1 < floodRainRange) {
                    Wu[i + 1] = 0.0;
                    //                        }

                    // 如果超过下层蒸散发能力
                    if (Wl[i] >= C * Lm) {
                        El[i] = (EPt - Eu[i]) * Wl[i] / Lm;
                        Ed[i] = 0.0;

                        /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                        //                            if (i + 1 < floodRainRange) {
                        Wl[i + 1] = Wl[i] - El[i];
                        Wd[i + 1] = Wd[i];
                        //                            }
                    }
                    // 如果没超过下层蒸散发能力，分两种情况讨论是否超过深层蒸散发能力，以C*(EPt-Eu[i])为分界
                    else {
                        // 如果大于，下层按能力蒸散发，深层无蒸散发
                        if (Wl[i] >= C * (EPt - Eu[i])) {
                            El[i] = C * (EPt - Eu[i]);
                            Ed[i] = 0.0;

                            /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                            //                                if (i + 1 < floodRainRange) {
                            Wl[i + 1] = Wl[i] - El[i];
                            Wd[i + 1] = Wd[i];
                            //                                }
                        }
                        // 如果小于，下层全部都蒸散发，深层有蒸散发
                        else {
                            El[i] = Wl[i];
                            Ed[i] = C * (EPt - Eu[i]) - Wl[i];

                            /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                            //                                if (i + 1 < floodRainRange) {
                            Wl[i + 1] = 0.0;
                            Wd[i + 1] = Wd[i] - Ed[i];
                            //                                }
                        }
                    }
                }
            }

            /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
            //                if (i + 1 < floodRainRange) {
            //得到下一时段上、下、深层张力水蓄量初值，限定阈值
            // 下一个时段初值不能小于零
            if (Wu[i + 1] < 0.0)
                Wu[i + 1] = 0.0;
            if (Wl[i + 1] < 0.0)
                Wl[i + 1] = 0.0;
            if (Wd[i + 1] < 0.0)
                Wd[i + 1] = 0.0;

            //            // 下一个时段初值不能大于最大张力水容量
            //            if (Wu[i + 1] > Um)
            //                Wu[i + 1] = Um;
            //            if (Wl[i + 1] > Lm)
            //                Wl[i + 1] = Lm;
            //            if (Wd[i + 1] > Math.max(0.0, Dm)) {
            //                Wd[i + 1] = Math.max(0.0, Dm);
            //            }
        }
    }

    /**
     * 分水源模型
     *
     * @return 地表水、壤中水、地下水
     */
    private void freeTank() {
        S = new double[floodRainRange + 1];// 存每一时段初产流面积上的自由水蓄量深
        FR = new double[floodRainRange + 1];// 存每一时段的产流面积（认为时段内的产流面积不变）
        //        S[0] = 0.0;// 假定初始自由水蓄量为0
        //        FR[0] = 0.001; // 假定初始产流面积为0.001
        S[0] = WInitialDay[3]; // 设定产流面积初始自由水蓄量
        FR[0] = WInitialDay[4]; // 设定初始产流面积系数

        // 参数KI和KG通常以日尺度定义，需要转换为时段尺度
        double D = 24 / timeInterval;
        double KID = (1 - Math.pow((1 - (KI + KG)), 1.0 / D)) / (1 + KG / KI);
        double KGD = KID * KG / KI;

        double EPt; // 蒸散发，等于Kc * E（实测蒸发，即ES[]或时段）
        double PEt; // PEt为净雨，等于降雨减去蒸发

        for (Integer i = 0; i < floodRainRange; i++) {

            // 蒸发计算
            if (floodE == null || floodE.length != floodRainRange) {
                EPt = Kc * calcEvaporationPotential(floodTm[i], timeInterval);
            } else {
                EPt = Kc * floodE[i];
            }
            //EPt = Kc * calcEvaporationPotential(floodTm[i], floodDrp[i], timeInterval);

            // PEt为降雨减去蒸发，可正可负
            PEt = floodDrp[i] - EPt;

            // 范围控制
            S[i] = S[i] < 0 ? 0 : S[i];

            // 有产流量才需要分段计算
            if (PEt > 0) {

                // 通过径流系数计算产流面积,产流面积等于流域总径流量/净雨
                FR[i] = R[i] / PEt;

                // 此处一定要控制产流面积的范围，因为计算出的产流面积有可能不在合理范围内
                if (FR[i] > 1.0)
                    // 注意产流面积不能比1大
                    FR[i] = 1.0;
                if (FR[i] < 0.0)
                    // 注意产流面积不能比0小
                    FR[i] = 0.0;

                if (i != 0) {
                    // 由于产流面积的变化，需将上个时段末的自由水蓄量换算成本时段初对应产流面积上的深度值）
                    S[i] = S[i] * FR[i - 1] / FR[i];
                }

                // 以5mm分段，划分成多段进行计算
                double q = R[i] / FR[i]; // 产流面积上净雨
                double N = Math.floor(q / 5 + 1);
                double pe = q / N;// 每段的净雨深

                // 时段长改变，参数也要相应变化(测试)
                /*double KIDN = (1 - Math.pow((1 - (KID + KGD)), 1.0 / N)) / (1 + KGD / KID);
                double KGDN = KIDN * KGD / KID;*/

                double SMM = SM * (1 + EX);// 流域单点最大自由水蓄水容量

                // 产流面积上的单点最大自由水蓄水容量
                double SMMF;

                if (EX <= 0.0001) {
                    SMMF = SMM;
                } else {
                    SMMF = SMM * (1.0 - Math.pow(1.0 - FR[i], 1.0 / EX));
                }

                double SMF = SMMF / (1.0 + EX); //产流面积上平均自由水蓄水容量

                double s = S[i]; // 该时段初产流面积上的自由水蓄量
                double rs = 0.0; // 该时段流域上的地表产流
                double ri = 0.0; // 该时段流域上的土壤产流
                double rg = 0.0; // 该时段流域上的地下产流
                for (Integer j = 0; j < N; j++) {
                    if (s >= SMF) { // 若自由水蓄水库本身已蓄满
                        rs = pe * FR[i];
                        ri = SMF * KID * FR[i];
                        rg = SMF * KGD * FR[i];

                       /* ri = SMF * KIDN * FR[i];
                        rg = SMF * KGDN * FR[i];*/
                        s = SMF - (ri + rg) / FR[i];
                    } else {
                        double AU = SMMF * (1 - Math.pow(1 - s / SMF, 1 / (1 + EX)));
                        if (pe + AU < SMMF) {
                            rs = FR[i] * (pe - SMF + s + SMF * Math.pow(1 - (pe + AU) / SMMF, 1 + EX));
                            double SMt = pe + s - rs / FR[i]; //时段产流面积上扣除地表产流后的自由水蓄量
                            ri = SMt * KID * FR[i];
                            rg = SMt * KGD * FR[i];

                            /*ri = SMt * KIDN * FR[i];
                            rg = SMt * KGDN * FR[i];*/
                            s += pe - (rs + ri + rg) / FR[i];
                        } else {
                            rs = (pe + s - SMF) * FR[i];
                            ri = SMF * KID * FR[i];
                            rg = SMF * KGD * FR[i];

                            /*ri = SMF * KIDN * FR[i];
                            rg = SMF * KGDN * FR[i];*/
                            s = SMF - (ri + rg) / FR[i];
                        }
                    }
                    RS[i] += rs;
                    RI[i] += ri;
                    RG[i] += rg;
                }

                /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                //                if(i + 1 < floodRainRange){
                S[i + 1] = s;
                //                }
            } else {
                // 无产流量，产流面积与上个时段一样
                if (i != 0) {
                    FR[i] = FR[i - 1];
                    // 此处一定要控制产流面积的范围，因为计算出的产流面积有可能不在合理范围内
                    if (FR[i] > 1.0)
                        // 注意产流面积不能比1大
                        FR[i] = 1.0;
                    if (FR[i] < 0.0)
                        // 注意产流面积不能比0小
                        FR[i] = 0.0;
                }
                RS[i] = 0.0;
                RI[i] = KID * S[i] * FR[i];
                RG[i] = KGD * S[i] * FR[i];

                /*需保留最后一个时段末的状态值，去掉时段初值限制条件*/
                //                if(i + 1 < floodRainRange){
                S[i + 1] = S[i] * (1 - KID - KGD);
                if (S[i + 1] < 0.0)
                    S[i + 1] = 0.0;
                if (S[i + 1] > SM)
                    S[i + 1] = SM;
                //                }
            }
            RS[i] = RS[i] * (1 - Im) + RIm[i]; // 地表径流量（mm）
            RI[i] = RI[i] * (1 - Im);           // 壤中径流量（mm）
            RG[i] = RG[i] * (1 - Im);           // 地下径流量（mm）
        }
    }

    /**
     * 坡面汇流
     * 线性水库地面径流、壤中流、地下径流
     *
     * @return 地表、壤中流、地下径流汇流
     */
    private void LinerReservoir() {
        double U = Area /
                (3.6 * timeInterval); //R 单位转换系数（km2/h = 1000000 m2/ 3600 s = 1000m2/3.6s，km2·mm/h = 1000000 * 0.001 m3/
        // 3600 s = m3/3.6s）

        // 参数CI和CG通常以日尺度定义，需要转换为时段尺度
        double D = 24 / timeInterval;
        double CID = Math.pow(CI, 1.0 / D);
        double CGD = Math.pow(CG, 1.0 / D);

        /*
        QS[0] = (RS[0] + RIm[0]) * (1 - CS) * U;
        QI[0] = RI[0] * (1 - CI) * U;
        QG[0] = RG[0] * (1 - CG) * U;
        */

        // 新增逻辑
        //        QS[0] = (RS[0] + RIm[0]) * U;
        QS[0] = RS[0] * U;
        QI[0] = qip + RI[0] * (1 - CID) * U;
        QG[0] = qgp + RG[0] * (1 - CGD) * U;
        QSS0[0] = qsp;
        QIS0[0] = qip;
        QGS0[0] = qgp;
        //        QS[0] = qsp;
        //        QI[0] = qip;
        //        QG[0] = qgp;

        for (Integer i = 1; i < floodRainRange; i++) {
            //            QS[i] = QS[i - 1] * CS + (RS[i] + RIm[i]) * (1 - CS) * U;
            //            QS[i] = (RS[i] + RIm[i]) * U;
            QS[i] = RS[i] * U;
            QI[i] = QI[i - 1] * CID + RI[i] * (1 - CID) * U;
            QG[i] = QG[i - 1] * CGD + RG[i] * (1 - CGD) * U;
            QSS0[i] = QS[i - 1];
            QIS0[i] = QI[i - 1];
            QGS0[i] = QG[i - 1];
        }

        QSS0[floodRainRange] = QS[floodRainRange - 1];
        QIS0[floodRainRange] = QI[floodRainRange - 1];
        QGS0[floodRainRange] = QG[floodRainRange - 1];

        //        for (Integer i = 1; i < floodRainRange; i++) {
        ////            QS[i] = QS[i - 1] * CS + (RS[i] + RIm[i]) * (1 - CS) * U;
        //            QS[i] = (RS[i] + RIm[i]) * U;
        //            QI[i] = QI[i - 1] * CID + RI[i] * (1 - CID) * U;
        //            QG[i] = QG[i - 1] * CGD + RG[i] * (1 - CGD) * U;
        //        }
    }

    /**
     * 河网汇流
     * 滞时法
     *
     * @param Q 单元面积河网总入流
     * @return 河网汇流-单元面积出口流量（m3/s）
     */
    private double[] TimeLag(double[] Q) {

        double[] Qf = new double[Q.length];
        Integer T = Integer.valueOf((int) L);
        if (T <= 0) {
            T = 0;
            for (Integer i = 0; i < Q.length; i++) {
                if (i == 0) {
                    //                    Qf[0] = (1 - CS) * Q[0];
                    Qf[i] = QSInitial[0];
                } else {
                    Qf[i] = CS * Qf[i - 1] + (1 - CS) * Q[i];
                }
            }
        } else {
            for (Integer i = 0; i < Q.length; i++) {
                /*
                if (i == 0) {
                    Qf[0] = 0.0;
                } else if (i < T) {
                    Qf[i] = CS * Qf[i - 1];
                    */
                if (i < T) {
                    if (i < Qf.length) {
                        if (i < QSInitial.length) {
                            Qf[i] = QSInitial[i];
                        } else {
                            Qf[i] = Qf[i - 1];
                        }
                    }
                } else {
                    Qf[i] = CS * Qf[i - 1] + (1 - CS) * Q[i - T];
                }
            }
        }

        // 新增逻辑:存储每一步的初值
        QS0 = new double[Q.length + 1][QSInitial.length];
        for (int i = 0; i < QSInitial.length; i++) {
            QS0[0][i] = QSInitial[i];
        }
        for (int i = 1; i <= Q.length; i++) {
            if (i < QSInitial.length) {
                for (int j = 0; j < QSInitial.length; j++) {
                    if (i + j < QSInitial.length) {
                        QS0[i][j] = QSInitial[j + i];
                    } else {
                        QS0[i][j] = Q[(i + j) - QSInitial.length];
                    }
                }
            } else {
                for (int j = 0; j < QSInitial.length; j++) {
                    QS0[i][j] = Q[i - QSInitial.length + j];
                }
            }
        }

        //        QS0 = new double[Q.length][QSInitial.length];
        //        for(Integer i = 0; i < Q.length; i++){
        //            for(Integer j = 0; j < QSInitial.length; j++){
        //                if(i + j < Q.length){
        //                    QS0[i][j] = Q[i + j];
        //                }else{
        //                    QS0[i][j] = QS0[i][j - 1];
        //                }
        //            }
        //        }

        return Qf;
    }

    /**
     * 河道演进
     * 分段马斯京根
     *
     * @param Qr 河网汇流-单元面积出口流量（m3/s）
     * @return 河道演进-单元面积出口流量（m3/s）
     */
    private double[] muskingum(double[] Qr) {
        //double[] QC = new double[Qr.length];
        double[] QC = new double[N];//河段下断面时段初始流量
        double[] QRouting = new double[Qr.length];//河道演进结果
        QX0 = new double[Qr.length + 1][N]; // 新增逻辑:存储每一步的初值
        for (Integer j = 0; j < N; j++) {
            if (j < QXInitial.length) {
                QC[j] = QXInitial[j];
            } else {
                QC[j] = QC[j - 1];
            }
            QX0[0][j] = QC[j];
        }

        double CD = 0.5 * timeInterval + KE - KE * XE;
        double C0 = (0.5 * timeInterval - KE * XE) / CD;
        double C1 = (0.5 * timeInterval + KE * XE) / CD;
        double C2 = 1 - C0 - C1;

        double QI1 = 0, QI2 = 0, QO1 = 0, QO2 = 0;//河段上、下断面时段初始流量和时段末流量
        if (C0 >= 0.0 && C2 >= 0.0) {// C0和C2需大于等于0，才满足马斯京根法的适用条件
            for (Integer i = 0; i < Qr.length; i++) {
                if (i == 0) {
                    QI1 = Qr[i];
                } else {
                    QI1 = Qr[i - 1];
                }
                QI2 = Qr[i];

                // 分段计算
                if (N > 0) {
                    for (Integer j = 0; j < N; j++) {
                        QX0[i][j] = QC[j];// 新增逻辑:存储每一步的初值
                        QO1 = QC[j];
                        QO2 = C0 * QI2 + C1 * QI1 + C2 * QO1;
                        QI1 = QO1;
                        QI2 = QO2;
                        QC[j] = QO2;
                    }
                } else {
                    // 新增逻辑:存储每一步的初值
                    QX0[i] = QXInitial;

                    QO2 = Qr[i];
                }

                /*
                if (i != 0) {
                    QC[i] = C0 * Qr[i] + C1 * Qr[i - 1] + C2 * QC[i - 1];
                } else {
                    QC[i] = Qr[i];
                }
                 */

                if (QO2 < 0.0001f)
                    QO2 = 0.0f;
                QRouting[i] = QO2;
            }

            for (Integer j = 0; j < N; j++) {
                QX0[Qr.length][j] = QC[j];
            }
        } else {// 否则，就不考虑河道演进了
            System.arraycopy(Qr, 0, QRouting, 0, Qr.length);
        }

        return QRouting;
    }

    private List<String> toStringList(double[] array0) {
        String[] array = new String[array0.length];
        for (int i = 0; i < array.length; i++) {
            array[i] = String.valueOf(array0[i]);
        }
        return Arrays.asList(array);
    }

    private List<String> toStringList(double[][] array0) {
        String[] array = new String[array0.length];
        for (int i = 0; i < array.length; i++) {
            array[i] = Arrays.toString(array0[i]);
        }
        return Arrays.asList(array);
    }

    @Override
    public void init() {

    }

    @Override
    public PredictResult predict() {
        PredictResult result = new PredictResult();
        Map<String, double[]> map = XAJmodel();

        double[][] RSim = new double[3][];
        //流域地表产流
        RSim[0] = map.get("RS");
        //流域土壤产流
        RSim[1] = map.get("RI");
        //流域地下产流
        RSim[2] = map.get("RG");

        //流域面积总产流，流域地表产流 RS[j] +  流域土壤产流RI[j] + 流域地下产流RG[j];
        result.setRunoffSim(map.get("runoffSim"));
        result.setRSim(RSim);
        //河道演进结果
        result.setQSim(map.get("QSim"));

        // 更新状态
        List<String> WUP = toStringList(Wu);
        List<String> WLP = toStringList(Wl);
        List<String> WDP = toStringList(Wd);
        List<String> SP = toStringList(S);
        List<String> FRP = toStringList(FR);
        List<String> QXSIG = toStringList(QX0);
        List<String> QSIG = toStringList(QS0);
        List<String> QSP = toStringList(QSS0);
        List<String> QIP = toStringList(QIS0);
        List<String> QGP = toStringList(QGS0);

        Map<String, List<String>> status = new HashMap<>();
        //上层张力水含量
        status.put("WUP", WUP);
        //下层张力水含量
        status.put("WLP", WLP);
        //深层张力水含量
        status.put("WDP", WDP);
        //产流面积自由水深
        status.put("SP", SP);
        //产流面积系数
        status.put("FRP", FRP);
        //河段下断面各时段初始流量
        status.put("QXSIG", QXSIG);
        //各时段初始总基流
        status.put("QSIG", QSIG);
        //地表坡面汇流
        status.put("QSP", QSP);
        //坡面汇流-壤中流
        status.put("QIP", QIP);
        //地下坡面汇流
        status.put("QGP", QGP);
        result.setNewStatus(status);
        return result;
    }

    //    /**
    //     * 土壤含水量初值计算
    //     *
    //     * @return 初始土壤含水量
    //     */
    //    private double[] calcAntecedentSoilMoisture() {
    //        double[] WFinalDay = new double[3];// 该方法的返回值W0，包括Wu0，Wl0，Wd0
    //        double Wmm = Wm * (1.0 + B) / (1.0 - Im); // 流域单点最大蓄水容量
    //        double[] WuDay = new double[preFloodRainRange + 1]; // 上层张力水蓄量，随时段变化
    //        double[] WlDay = new double[preFloodRainRange + 1]; // 下层张力水蓄量，随时段变化
    //        double[] WdDay = new double[preFloodRainRange + 1]; // 深层张力水蓄量，随时段变化
    //
    //        WuDay[0] = WInitialDay[0]; // 上层张力水蓄量初值
    //        WlDay[0] = WInitialDay[1]; // 下层张力水蓄量初值
    //        WdDay[0] = WInitialDay[2]; // 深层张力水蓄量初值
    //
    //        double[] EuDay = new double[preFloodRainRange]; // 上层蒸散发量，随时段变化
    //        double[] ElDay = new double[preFloodRainRange]; // 下层蒸散发量，随时段变化
    //        double[] EdDay = new double[preFloodRainRange]; // 深层蒸散发量，随时段变化
    //
    //        double[] Runoff = new double[preFloodRainRange];// 时段总产流
    //
    //        double EPt; // EPt为蒸散发能力，等于Kc * EM（水面蒸发，即dayEvapor[]）
    //        double PEt; // PEt为净雨，等于降雨减去蒸发
    //        double Wmt; // 某个时段初始的流域平均张力水蓄量
    //        double At; // 与Wmt值相对应的纵坐标值
    //        double WPt; //时段末上层张力水蓄量
    //        double WWPt;//时段末上层和中层张力水蓄量
    //
    //        if (Dm < 0) {
    //            System.out.println("Dm:" + Dm);
    //        } else {
    //            for (Integer i = 0; i < preFloodRainRange; i++) {
    //                Wmt = WuDay[i] + WlDay[i] + WdDay[i];
    //                if (Wmt < 0.0) {
    //                    System.out.println("W:" + Wmt);
    //                }
    //                if (Wmt > Wm) {
    //                    Wmt = Wm;
    //                }
    //                //蒸发能力计算
    //                EPt = Kc * calcEvaporationPotential(preFloodTm[i], preFloodDrp[i], preTimeInterval);
    //                PEt = preFloodDrp[i] - EPt; // PEt为降雨减去蒸发，可正可负
    //
    //                //如果有净雨产生，上、下、深层张力水蓄量主要受到降雨影响，是加法 *********/
    //                if (PEt >= 0.0) {
    //                    //计算透水面积产流
    //                    At = Wmm * (1.0 - Math.pow((1.0 - Wmt / Wm), 1.0 / (1.0 + B)));
    //                    if (PEt + At <= Wmm) {
    //                        Runoff[i] = PEt + Wmt - Wm + Wm * Math.pow((1 - (PEt + At) / Wmm), 1 + B);
    //                    } else {
    //                        Runoff[i] = PEt + Wmt - Wm;
    //                    }
    //
    //                    //各层蒸散发初值
    //                    EuDay[i] = EPt;
    //                    ElDay[i] = 0.0;
    //                    EdDay[i] = 0.0;
    //
    //                    /********* 分情况讨论上、下、深层张力水蓄量是否饱和，计算时段末WuDay、WlDay、WdDay *********/
    //                    // 如果上层未饱和
    //                    WPt = WuDay[i] + PEt - Runoff[i]; // PEt里面已经考虑蒸发了
    //                    if (WPt < Um) {
    //                        WuDay[i + 1] = WPt;
    //                        WlDay[i + 1] = WlDay[i];
    //                        WdDay[i + 1] = WdDay[i];
    //                    }
    //                    // 如果上层饱和，分情况讨论下层是否饱和
    //                    else {
    //                        // 在上层饱和前提下，如果下层未饱和
    //                        WWPt = WuDay[i] + WlDay[i] + PEt - Runoff[i];
    //                        if (WWPt < Um + Lm) {
    //                            WuDay[i + 1] = Um;
    //                            WlDay[i + 1] = WWPt - Um;
    //                            WdDay[i + 1] = WdDay[i];
    //                        }
    //                        // 在上层饱和前提下，如果下层饱和
    //                        else {
    //                            WuDay[i + 1] = Um;
    //                            WlDay[i + 1] = Lm;
    //                            WdDay[i + 1] = Wmt + PEt - Runoff[i] - WuDay[i + 1] - WlDay[i + 1];
    //                        }
    //                    }
    //                }
    //                //如果没净雨产生，上、下、深层张力水蓄量主要受到蒸发影响，是减法 *********/
    //                else if (PEt < 0.0) {
    //                    Runoff[i] = 0.0;
    //                    /********* 分情况讨论上、下、深层蒸散发是否达到能力，计算时段末WuDay、WlDay、WdDay *********/
    //                    // 如果超过蒸散发能力，即WuDay[i]+Rain[i]-EPt>=0，仅上层有蒸发
    //                    if (WuDay[i] + PEt >= 0.0) {
    //                        EuDay[i] = EPt;
    //                        ElDay[i] = 0.0;
    //                        EdDay[i] = 0.0;
    //
    //                        WuDay[i + 1] = WuDay[i] + PEt; // PEt=Rain[i]-EPt，这里PEt是负值
    //                        WlDay[i + 1] = WlDay[i];
    //                        WdDay[i + 1] = WdDay[i];
    //                    }
    //                    // 如果没超过蒸散发能力，即WuDay[i]+Rain[i]-EPt<0，分情况讨论是否超过下层蒸散发能力
    //                    else {
    //                        EuDay[i] = WuDay[i] + preFloodDrp[i];
    //                        WuDay[i + 1] = 0.0;
    //
    //                        // 如果超过下层蒸散发能力
    //                        if (WlDay[i] >= C * Lm) {
    //                            ElDay[i] = (EPt - EuDay[i]) * WlDay[i] / Lm;
    //                            EdDay[i] = 0.0;
    //
    //                            WlDay[i + 1] = WlDay[i] - ElDay[i];
    //                            WdDay[i + 1] = WdDay[i];
    //                        }
    //                        // 如果没超过下层蒸散发能力，分两种情况讨论是否超过深层蒸散发能力，以C*(EPt-EuDay[i])为分界
    //                        else {
    //                            // 如果大于，下层按能力蒸散发，深层无蒸散发
    //                            if (WlDay[i] >= C * (EPt - EuDay[i])) {
    //                                ElDay[i] = C * (EPt - EuDay[i]);
    //                                EdDay[i] = 0.0;
    //
    //                                WlDay[i + 1] = WlDay[i] - ElDay[i];
    //                                WdDay[i + 1] = WdDay[i];
    //                            }
    //                            // 如果小于，下层全部都蒸散发，深层有蒸散发
    //                            else {
    //                                ElDay[i] = WlDay[i];
    //                                EdDay[i] = C * (EPt - EuDay[i]) - WlDay[i];
    //
    //                                WlDay[i + 1] = 0.0;
    //                                WdDay[i + 1] = WdDay[i] - EdDay[i];
    //                            }
    //                        }
    //                    }
    //                }
    //                //得到下一时段上、下、深层张力水蓄量初值，限定阈值
    //                // 下一个时段初值不能小于零
    //                if (WuDay[i + 1] < 0.0)
    //                    WuDay[i + 1] = 0.0;
    //                if (WlDay[i + 1] < 0.0)
    //                    WlDay[i + 1] = 0.0;
    //                if (WdDay[i + 1] < 0.0)
    //                    WdDay[i + 1] = 0.0;
    //
    //                // 下一个时段初值不能大于最大张力水容量
    //                if (WuDay[i + 1] > Um)
    //                    WuDay[i + 1] = Um;
    //                if (WlDay[i + 1] > Lm)
    //                    WlDay[i + 1] = Lm;
    //                if (WdDay[i + 1] > Dm)
    //                    WdDay[i + 1] = Dm;
    //
    //                WFinalDay[0] = WuDay[i + 1];
    //                WFinalDay[1] = WlDay[i + 1];
    //                WFinalDay[2] = WdDay[i + 1];
    //            }
    //        }
    //
    //        return WFinalDay;
    //    }
}